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Further developments are introduced in the theory of finite size errors in quantum many-body simulations of 
extended systems using periodic boundary conditions. We show that our recently introduced Model Periodic 
Coulomb interaction [A. J. Williamson et al., Phys. Rev. B 55, R4851 (1997)] can be applied consistently to all 
Coulomb interactions in the system. The Model Periodic Coulomb interaction greatly reduces the finite size errors 
in quantum many-body simulations. We illustrate the practical application of our techniques with Hartree-Fock 
and variational and diffusion quantum Monte Carlo calculations for ground and excited state calculations. We 
demonstrate that the finite size effects in electron promotion and electron addition/subtraction excitation energy 
calculations are very similar. 

PACS: 71.10.-w, 71.15.-m, 71.15.Nc 



I. INTRODUCTION 

Most simulations of extended quantum systems are 
performed using finite simulation cells. This introduces 
"finite size errors" which are one of the major problems 
limiting the application of accurate many-body tech- 
niques to extended systems. The standard method of re- 
ducing the finite size errors is to apply periodic boundary 
conditions, but important finite size errors often remain. 
In this paper we present further developments of the the- 
ory of finite size effects in quantum many-body simula- 
tions subject to periodic boundary conditions. Our mo- 
tivation is to understand and reduce the finite size ef- 
fects encountered in quantum Monte Carlo simulations, 
although the techniques described here are of wide gen- 
erality and can be readily applied to other many-body 
electronic structure methods. 

Quantum Monte Carlo (QMC) methods in the varia- 
tional [Q and diffusion |^],[| forms are capable of yielding 
highly accurate results for correlated electron systems. 
These methods are very promising because (i) electron 
correlations are included essentially without approxima- 
tion and (ii) the methods scale favourably with system 
size. Nevertheless, for realistic systems the cost of these 
calculations remains large, and it would be highly de- 
sirable to reduce the finite size errors so that accurate 
results can be obtained using small simulation cells. We 
stress that all quantum many-body calculations which 
use periodic boundary conditions to model extended sys- 
tems suffer from finite size effects and that the ideas dis- 
cussed in this paper are relevant whenever long ranged 
interactions are involved. 

Finite size errors in many-body calculations have tra- 
ditionally been corrected for by extrapolation techniques 
and/or by using the results of more approximate calcu- 
lations, such as those based on density functional theory 
(DFT). While extrapolation techniques can certainly be 
successful they are very costly. In addition, the finite 



size errors in periodic boundary condition DFT calcula- 
tions are significantly different from those in the standard 
formulations of quantum many-body techniques such as 
QMC, and therefore a simple application of a DFT finite 
size correction may not lead to accurate results. We have 
understood the reason for this difference and have found 
a way to reduce the finite size errors in quantum many- 
body simulations using a new "Model Periodic Coulomb" 
(MPC) interaction, which has the additional advantage 
that the residual finite size effects are reasonably well 
described by standard DFT calculations. The accuracy 
can be further increased by using an extrapolation proce- 
dure, but the extrapolation corrections are considerably 
reduced and can therefore be evaluated using a smaller 
range of system sizes. 

The layout of this paper is as follows. In section || we 
describe the Hamiltonian within periodic boundary con- 
ditions, while in section 111 we discuss various finite size 
correction and extrapolation procedures. In section IV 



we review the "independent particle finite size effects" , 
showing how our k-space sampling techniques are related 
to those used in mean-field theories. In section |v| we 
introduce our MPC interaction for reducing finite size 
effects in periodic systems and show that it can be ap- 
plied to all the Coulomb interactions. We presen t test s 
of the MPC in teract ion within HF theor y (sec tion VIA), 
VMC (section [vf|), and DMC (section |VIC| ). The lat- 
ter section includes DMC results for a system with 1000 
electrons, which is the largest number in any DMC cal- 
culation to date. In section VII we discuss finite size er- 
rors present in calculations of excitation energies. Tests 



within HF theory are presented in section VII A, while 
in section VII B we present VMC results for the "op- 
tical absorption" and "photoemission" gaps, the latter 
being the first such calculations for a three-dimensional 
periodic system. Calculations with up to 512 electrons 
demonstrate that the finite size errors in the "optical ab- 
sorption" and "photoemission" gaps are similar and that 



1 



the finite size effects are quite small even for a 64-electron 
simulation cell. We draw our conclusions in section VIII. 



III. FINITE SIZE CORRECTION AND 
EXTRAPOLATION FORMULAE 



II. THE HAMILTONIAN WITHIN PERIODIC 
BOUNDARY CONDITIONS 

The many-body Hamiltonian for a system of electrons 
at positions and static ions at positions x Q is Q] 

H = E-^Vf + ^5> Q (r 4 ,x Q ) 

i i a 

+ lj2J2 V ( ri ' r i"> + lYl E ^(XoX/s) ■ (1) 
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An infinite system is normally modeled by a finite simu- 
lation cell subject to periodic boundary conditions. The 
model interaction terms, v a , v, and v a p, are chosen such 
that the potential energy of the model system, which 
involves only the positions of the particles in the finite 
simulation cell, mimics the potential energy of the infi- 
nite system as closely as possible. Since the potential 
energy of the infinite system depends on the positions 
of all the charges in the solid, only a few of which are 
included in the simulation, the model interaction energy 
is approximate even in crystalline solids. To enforce the 
periodic boundary conditions the functions v a , v, and 
v a /3 must be invariant under the translation of either ar- 
gument by a member of the set of translation vectors 
of the simulation cell lattice, {R}. The standard ap- 
proach is to choose the model Hamiltonian such that the 
full potential energy of Eq. [j], evaluated by summing the 
model interactions between all pairs of particles in the 
simulation cell, equals the potential energy per cell of an 
infinite array of identical copies of the simulation cell. 
However, even when we restrict each simulation cell in 
the array of copies to be overall charge neutral, the sum 
of inter -particle Coulomb 1/r interactions is only condi- 
tionally convergent, H and to define this model interac- 
tion uniquely the boundary conditions at infinity must 
be specified. The standard procedure is to define the po- 
tential by solving Poisson's equation subject to periodic 
boundary conditions, in which case the model interaction 
is the Ewald interaction. || For two electrons separated 
by r, the Ewald interaction is 



E 

{G}^0 



4-7T 

^2 exp [-G 2 /4k 2 + zG.r] 



7T ^ er f C ( K l r 

{R} 1 



Rl 



Rl 



(2) 



where fl is the volume of the simulation cell, {G} is the 
set of reciprocal space translation vectors of the simula- 
tion cell lattice, and k is a positive but otherwise arbi- 
trary constant. 



The basic idea behind finite size correction formulae is 
to write the energy for the infinite system as 



En + (Eoo — En) 



(3) 



where the subscript denotes the system size. A highly 
accurate many-body calculation for the ./V-particle sys- 
tem is then performed to obtain an approximation for 
En, and the correction term in brackets is approximated 
using a much less expensive scheme which can be applied 
to very large systems. 

In an extrapolation procedure the energy En is calcu- 
lated for a range of system sizes and fitted to a chosen 
functional form which contains parameters. Correction 
and extrapolation procedures can be combined to give an 
expression 



E N + - E' N ) + F(N) , 



(4) 



where the prime indicates that a less expensive scheme 
is used and F(N) is an extrapolation function. Clearly, 
the optimal form of F(N) depends on the method used 
to calculate the correction term. The practical difference 
between correction and extrapolation is that correction 
requires a single calculation of En using the accurate 
(and normally very costly) many-body technique, while 
extrapolation requires several such calculations for dif- 
ferent values of N and a subsequent fit to the chosen 
functional form F(N). The extrapolation procedure is 
costly because it involves more calculations and is prone 
to inaccuracy because one has to perform a fit with only a 
few data points. In designing a correction/extrapolation 
procedure one therefore tries to make the extrapolation 
term as small as possible. 

Candidate methods for evaluating the correction term 
include Kohn-Sham DFT and HF theory. The most 
convenient methods arc the local density approximation 
(LDA) to DFT, or extensions such as the Generalized 
Gradient Approximation. These methods are very widely 
applied in periodic boundary conditions calculations and 
are computationally inexpensive, while retaining a real- 
istic description of the system. Within an independent 
particle theory, such as Kohn-Sham DFT, calculations 
for periodic systems are normally performed by solving 
within the primitive unit cell. To obtain the correct 
result for the infinite system it is necessary to integrate 
over k-space, and the integral is normally approximated 
by a sum over a finite set of k-points. A determinant 
formed from the occupied orbitals at a single k-point in 
the first Brillouin zone (BZ) of the primitive unit cell is 
a many-body wave function for a simulation cell of the 
size of the primitive unit cell. Adding a second k-point 
doubles the size of the determinant and is equivalent to 
doubling the size of the simulation cell, etc. 



2 



An important subtlety arises when finite size correc- 
tions derived from an independent particle theory such 
as LDA-DFT are used to correct the results of a true 
many-body method such as QMC. Suppose the many- 
body calculation is performed using the Ewald interac- 
tion, so that the Coulomb energy is that of an infinite 
periodic array of copies of the simulation cell. Given that 
the solid we are trying to model is crystalline, and hence 
the charge density is truly periodic, the Ewald interac- 
tion gives a good description of the classical Coulomb 
or Hartree energy. However, because the electronic po- 
sitions are mirrored exactly in every copy of the simula- 
tion cell, the electronic correlations are also forced to be 
periodic, and the exchange-correlation (XC) energy cor- 
responds to a system with a periodic XC hole. This un- 
physical approximation is particularly inaccurate when 
the simulation cell is small. In Kohn-Sham DFT cal- 
culations the XC energy is evaluated using a standard 
functional such as the LDA of Perdew and Zunger ||, 
which was obtained by fitting to the results of DMC cal- 
culations for jellium. |1 There is an important difference, 
however, in that the DMC results were extrapolated to 
the infinite-system-size limit before the fit was made. 

The consequences of this difference are nicely illus- 
trated by considering many-body and LDA calculations 
for jellium, in which the charge density is uniform. The 
LDA gives the XC energy for the infinite system irrespec- 
tive of the size of the simulation cell, but the XC energy 
obtained from a many-body simulation using the Hamil- 
tonian of Eq. [I] with the Ewald interaction gives an XC 
energy which depends on the size and shape of the simu- 
lation cell. Consequently, evaluating the correction term 
in Eq. ^ within the LDA does not give a good approx- 
imation to the finite size correction in the many-body 
simulation and a significant extrapolation term remains. 

The issue of finite size corrections to both the kinetic 
and interaction energies has been addressed by Ceperley 
and coworkers. |Ic|-[T3]| Their approach is to add separate 
extrapolation terms for the kinetic and interaction en- 
ergies. In their work on hydrogen solids, Ceperley and 
Alder performed DMC calculations for a number of 
different system sizes and fitted to the formula 

E UMC ^ E DMC + a{ToQ Tn) + b_ (5) 

where a and b are parameters, and T is the kinetic en- 
ergy of the non-interacting electron gas. The b/N term 
accounts for the finite size effects arising from the in- 
teraction energy and the difference of the parameter a 
from unity accounts for the difference between the kinetic 
energies of the interacting and non-interacting systems. 
Normally a > 1 because the interacting kinetic energy 
is larger than the non-interacting kinetic energy. Engel 
et al. ]l4| used the following formula for inhomogeneous 
systems: 



- + a(E^ - E]f A ) + ^ , (6) 

which reduces to Eq. |^ for a homogeneous system. 

As mentioned above, a large part of the finite size er- 
ror in the interaction energy arises from the use of the 
Ewald interaction. [ p3||l6| ] Our approach to the problem 
of the finite size errors differs from that of Ceperley et 
al. in an important respect. We try to reduce the finite 
size effects within the many body calculation by modify- 
ing the interaction terms in the many-body Hamiltonian. 
The Ewald interaction is a periodic function which dif- 
fers from 1 jr in such a way that the sum of interactions 
between pairs of particles within one cell gives the exact 
Coulomb energy per cell of a periodic lattice of identi- 
cal cells. This ensures that the Ewald Hamiltonian gives 
the correct Hartree energy, but the deviations from 1/r 
give rise to a spurious contribution to the XC energy, 
corresponding to the periodic repetition of the XC hole 
discussed earlier. ]l5],[l6| Our solution is to change the 
many-body Hamiltonian so that the interaction with the 
XC hole is exactly 1/r (see section |v|) for any size and 
shape of simulation cell, without altering the form of the 
Hartree energy. This source of finite size error is therefore 
eliminated and Eq. with the correction term calculated 
within the LDA, gives a much better description of the 
remaining finite size errors. Greater accuracy can be ob- 
tained by adding an extrapolation term, but this term is 
much smaller than if the Ewald interaction is used. 

The idea of changing the Hamiltonian to reduce the 
finite size errors may seem strange at first, so it is worth 
discussing a simple example. Imagine that a periodic 
boundary condition QMC technique is being used to 
study an isolated molecule. If the molecule is placed 
at the center of the periodically repeated simulation 
cell, the calculated energy is that of an array of iden- 
tical molecules, including unwanted inter-molecular in- 
teractions. The results improve as the simulation cell 
is made larger, but the convergence is slow, especially 
for molecules with permanent dipole moments. A better 
solution is to cut off all Coulomb interactions between 
charges on different molecules, i.e., to replace the Ewald 
interactions by truncated Coulomb 1/r interactions act- 
ing only within the simulation cell. As long as the molec- 
ular wave function has decayed to zero before the sim- 
ulation cell boundary is reached, this procedure should 
give essentially exact results. The changes to the interac- 
tion to be discussed in this paper are a generalization of 
this approach to make it useful in simulating genuinely 
periodic systems such as crystals. In this case the very 
"strict" notion of periodicity built into the Ewald inter- 
action produces an artificial periodic replication of the 
XC hole as well as the required periodic replication of 
the charge density. Our modified Hamiltonian removes 
the effects of the unwanted extra periodicity just as the 
simple truncation removed the unwanted periodicity in 
the molecular example. 
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An alternative procedure is to evaluate the correction 
term in Eq. |^ using HF data. HF theory is an approx- 
imate method for solving the many-body Hamiltonian, 
and if we use the Ewald formula for the electron-electron 
interaction terms in both the many-body and HF theo- 
ries, the finite size error in the HF exchange energy will 
tend to cancel the finite size error in the many-body in- 
teraction energy. At first sight this appears to be an 
excellent solution to the finite size problem. However, 
this procedure gives too large a correction, presumably 
because the HF exchange hole is significantly different 
from the screened XC hole of the many-body system. 



IV. INDEPENDENT PARTICLE FINITE SIZE 
EFFECTS 



have described the theory of the "special k-points" for 
many-body calculations in our earlier work p9 20 and 



We call the correction term (E^ )A - E}f A ) the "in- 
dependent particle finite size effect" (IPFSE). The finite 
size effect arising from the particular model interaction 
used is called the "Coulomb finite size effect" (CFSE). 
In this section we discuss the IPFSE in more detail. The 
"finite size error" in an LDA calculation for a perfect 
crystal arises from errors in the BZ integration. In an 
LDA calculation the computational cost is proportional 
to the number of k-points used. However, in a QMC cal- 
culation the volume of the simulation cell is proportional 
to the number of k-points in the primitive BZ and the 
computational cost increases approximately as the cube 
of the volume of the simulation cell. This means that it 
is even more important to choose the k-points carefully 
in a QMC calculation so as to make the IPFSE as small 
as possible. 

Additional errors arise when using the supercell ap- 
proximation for systems which break translational sym- 
metry. Consider a crystal containing a point defect. In 
the supercell approach a finite simulation cell containing 
a single defect is repeated throughout space to form an 
infinite three-dimensional array. The supercell must be 
sufficiently large that the interaction between defects in 
different cells is negligible. This type of finite size effect 
is distinct from the BZ integration error because it per- 
sists even when the BZ integration is performed exactly. 
However, the correction to the defect formation energy 
from this type of finite size error should be described rea- 
sonably well by a (£^ DA — Ej^ A ) correction term, where 
E}^ A and -E^ DA are, respectively, the defect formation 
energies for a large simulation cell and a smaller one, each 
containing a single defect. 

A determinantal Bloch wave function suitable for use 
in a QMC calculation may be formed from a set of sin- 
gle particle orbitals at a single k-point, k s , in the Bril- 
louin zone of the simulation cell reciprocal lattice. The 
IPFSE can be greatly reduced for insulating systems by 
a careful choice of k s using the method of "special k- 
points" borrowed from band structure theory. PJ,[18| We 



here we concentrate on the practicalities of choosing k s . 
Baldereschi |l7j defined the "mean-value point", which 
is a k-point at which smooth periodic functions of wave 
vector accurately approximate their averages over the 
BZ, and clearly the Baldereschi mean-value point is a 
strong candidate for k s . However, we find it convenient 
to choose k s to be equal to half a translation vector of the 
simulation cell reciprocal lattice (k s = G/2), which al- 
lows the construction of real single-particle orbitals and 
hence is computationally more efficient. Some freedom 
is still left in the selection of k s , and we choose from the 
G/2 according to the symmetrized plane wave test of BZ 
integration quality introduced by Baldereschi. jlTj For 
example, for an fee simulation cell the half-reciprocal- 
lattice vectors correspond to the T, X and L points of 
the BZ of the simulation cell reciprocal lattice. For a 
crystal with the full cubic symmetry k s = L gives the 
best BZ integration and k s = T the worst. 

It is illuminating to relate these k-point schemes to 
the multi-point schemes used in LDA and HF calcula- 
tions. In LDA and HF calculations one normally sam- 
ples the BZ of the primitive unit cell, whereas in a QMC 
calculation one samples the BZ of the simulation cell. 
Suppose that the simulation cell has translation vectors 
Ni&i, N2&2, and N^a^, where the Ni are integers and the 
aj are vectors defining the primitive unit cell. When un- 
folded into the BZ of the primitive unit cell the k-point 
k s maps onto the regular mesh 



' 1 m 1 n > 

k; m „ = — bi + — b 2 + — b 3 + k s 

ivi iV2 iV3 



(7) 



where; = 0, . . .N t -1, m = 0, . . . N 2 -l, n = 0, . . . iV 3 -l, 
and the b^ are reciprocal to the a^. This mesh is of the 
same type as defined in the widely-used Monkhorst-Pack 
scheme for BZ integration, differing only by an off- 
set from the origin. The multi-point generalization of 
the Baldereschi scheme which can be used in LDA and 
HF calculations is now obvious: one chooses the offset 
k s to be equal to the mean-value point |T^] of the lat- 
tice defined by the translation vectors A^a,. As men- 
tioned above, for our purposes we prefer to choose k s to 
be equal to half a translation vector of the simulation 
cell reciprocal lattice. This choice gives a smooth and 
rapid convergence of the BZ integration with increasing 
values of the Ni p]( ]. For example, as mentioned above, 
for fee crystals we choose k s to be an L-point of the 
BZ of the lattice defined by the A^a.;. In this case the 
Monkhorst-Pack Jllj definition of the offset corresponds 
to taking an L-point for Ni even, but the T-point for N 
odd. The latter choice gives poor results and this prob- 
lem has prompted some researchers to avoid Monkhorst- 
Pack meshes with odd values of the N, but in fact odd 
values can be very efficient if the mesh is offset according 
to our prescription. 
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V. COULOMB FINITE SIZE EFFECTS 



where p is the electronic charge density and 



In this section we summarize the theory of the MPC 
interaction. More details and background are given in 
Refs. |H1 and @. In section [II we described how the 



(10) 



CFSE arises from the XC energy and the dependence of 
the Ewald interaction, ve , on the size and shape of the 
simulation cell. Expanding ve around zero separation 

gives capi 



, . 1 277 rp 

ve(t) = constant H 1 — — r D r 

r 30 



O 



, (8) 



where f2 is the volume of the simulation cell, and the ten- 
sor D depends on the shape of the simulation cell. (For a 
cubic cell D=I.) The constant term arises from the con- 
dition that the average of ve over the simulation cell is 
zero. The term of order r 2 and the higher order devia- 
tions from 1/r make the Ewald interaction periodic and 
ensure that the sum of interactions between the particles 
in the simulation cell gives the potential energy per cell 
of an infinite periodic lattice. These terms are responsi- 
ble for the spurious contribution to the XC energy, which 
is the source of the large finite size effect in many-body 
calculations using the Ewald interaction. [jl5],[l6) For cu- 
bic cells the interaction at short distances is larger than 
1/r and therefore the XC energy is more negative than 
it should be, and because the leading order correction 
is proportional to the inverse of the simulation cell vol- 
ume the error per electron is inversely proportional to 
the number of electrons in the cell. 

Clearly it is desirable to remove this spurious contri- 
bution to the XC energy, but we must remember that the 
Hartree energy is correctly evaluated using the Ewald in- 
teraction. The key requirements for a model Coulomb 
interaction giving small CFSEs in simulations with pe- 
riodic boundary conditions are therefore: (i) it should 
give the Ewald interaction for the Hartree terms, and 
(ii) it should be exactly 1/r for the interaction with the 
XC hole. Unfortunately, the only periodic solution of 
Poisson's equation for a periodic array of charges is the 
Ewald interaction, which obeys criterion (ii) only in the 
limit of an infinitely large simulation cell. We therefore 
abandon the use of Poisson's equation and the Ewald in- 
teraction, and instead use a "Model Periodic Coulomb" 
(MPC) interaction which satifies both criteria. This may 
seem a drastic step, but we point out that we are trying 
to model an infinite system by a finite one and that the 
inter-particle interaction used in the finite simulation cell 
must model the effects of all the charges in the infinite 
system. 

Our MPC interaction [fill can be written as 



+ J2 [ [M'i - r ) - /(r* ~ r )] P( r ) dr . ( 9 ) 

, JWS 



The definition of the cutoff Coulomb function, /, involves 
a minimum image convention whereby the inter-electron 
distance, r, is reduced into the Wigner-Seitz (WS) cell 
of the simulation cell lattice by removal of simulation cell 
lattice translation vectors, leaving a vector r m . This en- 
sures that -He-e has the correct translational and point 
group symmetry. The first term in Eq. |^ is a direct 
Coulomb interaction between electrons within the sim- 
ulation cell and the second term is a sum of potentials 
due to electrons "outside the simulation cell" . Note that 
the second term is a one-body potential similar to the 
Hartree potential. It depends on the electronic charge 
density, p, but is not a function of the separation of the 
electrons. 

The electron-electron contribution to the total energy 
is evaluated as the expectation value of H e — e with the 
many-electron wave function, <f> 1 minus a double counting 
term: 



- ~ / p{v)p{v') Mr - r') - /(r - r')] dr dr' . (11) 
1 Jws 

Evaluating the expectation value gives 

£ c -c = \ f p(r)p(r> E (r - r') dr dr' 
z Jws 

+ {[ H 2 E/( r ^-) n * dr * 

\Jws t>j 

- \ [ p(r)p(r')/(r-r')drdr') , (12) 
L Jws J 



'WS 

where the first term on the right hand side is the Hartree 
energy and the term in brackets is the XC energy. We 
can see immediately that the Hartree energy is calcu- 
lated with the Ewald interaction while the XC energy 
(expressed as the difference between a full Coulomb term 
and a Hartree term) is calculated with the cutoff inter- 
action /. 

The charge density p appears in Eqs. ||, |ll|, and [l2|. In 
QMC methods, the charge density is known with greatest 
statistical accuracy at the end of the calculation. This is 
not a serious complication for VMC simulations as the in- 
teraction energy may be evaluated at the end of the sim- 
ulation using the accumulated charge density. In DMC 
this is not possible because the local energy [FJ,|| is re- 
quired at eve ry step. We investigate this point further in 
section VI C . 

The CFSE may be viewed in another way [^5| . Con- 
sider a large cluster of identical simulation cells. Almost 
all possible configurations of the electrons within the 
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simulation cell have a net dipole moment. The dipole 
moments in the periodic replicas of the simulation cell 
are aligned, resulting in a depolarization field across the 
sample. (In real non-ferroelectric solids the dipoles from 
different regions are not aligned and the depolarization 
field is greatly reduced.) To make the electrostatic po- 
tential periodic deep within the cluster we need to apply 
a cancelling external electric field. In the limit of a very 
large cluster this is equivalent to using the Ewald inter- 
action. In a cubic system the total energy of interaction 
between the dipoles, which includes the interaction with 
the depolarization field, does not contribute to the 0(r 2 ) 
term in Eq. || |21|. For a cubic system the 0(r 2 ) term 
in Eq. |§| therefore arises solely from the interaction be- 
tween the dipoles and the external field. In a non-cubic 
system the 0(r 2 ) term contains contributions from both 
the dipole-dipole interactions and the interactions of the 
dipoles with the external field. The advantage of our 
MPC interaction is that it prevents the 0(r 2 ) terms con- 
tributing to the XC energy. 

If the simulation cell is not big enough to accommodate 
the true XC hole then the XC hole is squeezed into the 
cell and a finite size error ensues. This appears to be 
the largest finite size error in the interaction energy after 
the CFSE has been removed. Extrapolation currently 
appears to be the best method of correcting this error. 
When the XC hole is squeezed by the finite simulation cell 
the XC energy calculated using the MPC interaction is 
too negative. The effect of the 0(r 2 ) term in the Ewald 
interaction is to make it even more negative, however, 
and so the MPC interaction is still better than the Ewald 
interaction. 



A. Systems of electrons and nuclei 

In the previous sections we considered the electron- 
electron interaction only, and did not discuss the 
electron-nucleus and nucleus-nucleus interactions. If the 
MPC interaction is physically reasonable we should, how- 
ever, be able to apply it to all the interactions in the 
problem, not just the electron-electron terms. Here we 
apply the MPC interaction to a system of electrons and 
nuclei, showing that under certain common conditions 
the expressions simplify so that the electron-nucleus and 
nucleus-nucleus terms reduce to the Ewald form. 

Consider a simulation cell with periodic boundary con- 
ditions containing N electrons at positions and M nu- 
clei of charge Z a at positions x a . The wave function of 
the electrons and nuclei is \&({rj}, {x Q }), and the total 
charge density (electrons and nuclei) is px (r) ■ The inter- 
action energy calculated with the MPC interaction is 



Eint — - 



\ [ p T (r)p T (r')K(r-r')-/(r-r')]drdr' 
L Jws 
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n fe dr k II 7 dx 7 



(13) 



We now employ the adiabatic approximation to sepa- 
rate the electronic and nuclear dynamical variables: 



9({ ri }, {x Q }) = 0({ r< }; {x Q })0>({x Q }) , 



(14) 



where the {x a } appear only as parameters in <f>. To 
make further progress we must assume a form for the 
nuclear part of the wave function, The simplest 

assumption is that $ can be written as an appropri- 
ately symmetrized product of very strongly localized 
non-overlapping single-nucleus functions. Regardless of 
whether the product is antisymmetrized, symmetrized, 
or not symmetrized, Eq. [l^ reduces to 

Eint = \ ! p(r)p(r> E (r - r') - /(r - r')] drdr' 
1 Jws 

Jws l>3 

- / |0| 2 y2 y^ZaV-E^i -x Q )n fe dr fe 
Jws 



a>(3 



(15) 



where the x Q denote the centers of the single-nucleus 
functions and p is the electron density. Note that the 
electron-nucleus and nucleus-nucleus terms involve only 
the Ewald interaction and that the first two terms of 
Eq. |l5] correspond precisely to the electron-electron in- 
teraction of Eq. [l2| This result justifies the use of Eq. ^ 
for the electron-electron interactions while retaining the 
Ewald interaction for the electron-nucleus and nucleus- 
nucleus terms. 



VI. TESTS OF THE MPC INTERACTION 

A. Application to HF calculations 

We have tested the MPC interaction by performing a 
series of calculations on diamond structure silicon using 
fee simulation cells whose translation vectors are n times 
those of the primitive unit cell. In a previous publication 
we gave a few results from such tests [fl6|j , but here we 
present new tests and subject them to a more detailed 
analysis. 

In our first set of tests we compare LDA and HF re- 
sults. These tests are inexpensive, which allows us to 
study very large systems, and they are not subject to 
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statistical errors because Monte Carlo techniques are not 
involved. We consider simulation cells with n = 1, 2, 
3, 4, and 5, which contain 2, 16, 54, 128, and 250 ions, 
respectively. The Si 4+ ions were represented by norm- 
conserving non-local LDA pseudopotentials and the cal- 
culations were performed using a plane-wave basis set 
with a cutoff energy of 15 Ry. To facilitate comparison 
we evaluate the HF energy with the LDA orbitals, so that 
the energy differences arise solely from the difference be- 
tween the LDA XC energy and the HF exchange energy. 
We performed calculations using k s = L and k s = T. 
The HF energy evaluated with the MPC interaction is 

SHF = E4 / ^(r)V 2 <Mr)dr 
+ \ y '°( r )' ( r/ ) w E(i' - r') drdr' 

N 

- ^W^(r')/(r-r')0*(r')^(r)drdr' 



V ext (r)p(r) dr + £i_i 



(16) 



where Ei—i is the ion-ion energy calculated with the 
Ewald interaction. Note that the Hartree energy is eval- 
uated with the Ewald interaction while the exchange en- 
ergy is evaluated with the cutoff Coulomb interaction. 

In Fig. [I] we show the deviations of the LDA and HF 
energies from the fully converged values as a function of 
system size for k s = T wave functions. The LDA en- 
ergy converges smoothly with system size but for small 
system sizes the IPFSE error is large because of the T- 
point sampling. We show HF energies calculated with the 
Ewald and MPC interactions, with and without incor- 
porating the IPFSE corrections obtained from the LDA 
data. The data incorporating the IPFSE (filled symbols) 
show the residual CFSE errors. The IPFSE is positive 
while the CFSE is negative, in accord with the analysis of 
the CFSE given in section [v|. The IPFSE corrected data 
show that the CFSE for the MPC interaction is roughly 
half that for the Ewald interaction. 

In Fig. we show similar data for k s = L. The LDA 
energy converges rapidly and smoothly with system size, 
and therefore the IPFSE is small, which demonstrates the 
efficacy of the "special k-points" method. The IPFSE 
and CFSE errors tend to cancel and for k s = L sampling 
the HF data converge more rapidly without the IPFSE 
corrections. For n — 3, corresponding to a 54 atom sim- 
ulation cell, the LDA finite size error (IPFSE only) is 
0.011 eV per atom, which is very much smaller than the 
HF (Ewald) finite size error of -0.211 eV per atom or the 
equivalent HF (MPC) error of -0.071 eV per atom. As 
in the case of F-point sampling, the CFSE for the MPC 
interaction is roughly half that for the Ewald interaction. 

After applying the IPFSE correction obtained from the 
LDA data the HF results for k s = T and k s = L sampling 



are very similar. The correspondence of the IPFSE cor- 
rected data for the L— and T-points demonstrates that to 
a very good approximation the IPFSE and CFSE are in- 
dependent. Estimating the energy of the infinite system 
by averaging the energy over a set of k s vectors removes 
the IPFSE but does not remove the CFSE. 

We have fitted the CFSE errors from the filled data 
points in Figs. |l] and [| to the form b/N x , where b and x 
are parameters and N = 8n 3 is the number of electrons 
in the simulation cell. The fits give values of x in the 
region of unity for both the Ewald and MPC interactions. 
This extrapolation function is therefore suitable for both 
interactions, although the size of the extrapolation term 
is smaller for the MPC interaction. 

The Ewald and MPC interactions differ by an amount 
inversely proportional to the volume of the simulation 
cell. This means that although the energies per particle 
converge to the same value as the volume of the sim- 
ulation cell increases, the difference between the Ewald 
and MPC energies of the whole simulation cell tends to 
a finite constant as the simulation cell volume goes to 
infinity. We have evaluated this constant value for the 
fixed-LDA-orbital HF calculations described in this sec- 
tion by extrapolating the energy difference between the 
Ewald and MPC energies. This gives a value of approxi- 
mately 14 eV. 



B. Application to VMC 

In the VMC method §,§ the energy is calculated as 
the expectation value of the Hamiltonian, H , with a trial 
wave function, </>t, yielding a rigorous upper bound to the 
exact ground state energy. The Metropolis algorithm is 
used to generate electron configurations distributed ac- 
cording to 4>\i and the energy calculation is performed 
by averaging the local energy, (f)^ 1 Hcp^, over this distri- 
bution. 

Our trial wave functions are of the standard Slater- 
Jastrow type: 



b T = D^D 1 exp 



JY 



N 



(17) 



where there are N electrons in the simulation cell, x is 
a one-body function, u is a two-body correlation factor 
which depends on the relative spins of the two electrons, 
and and D^- are Slater determinants of up- and down- 
spin single-particle orbitals. The u functions were of the 
type described in |22|], while for the \ functions we used 
spherically symmetric functions centered on each atom. 
These x functions give significantly better results than 
the truncated Fourier series representation used in our 
earlier work. [^2| , [l6| The trial wave functions contained 
32 variable parameters, whose optimal values were ob- 
tained by minimizing the variance of the energy using 
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10,000-20,000 statistically independent electron configu- 
rations, which were regenerated several times during the 
minimization procedure. |^,^3| We used the same pseu- 
dopotential as in our LDA and HF calculations, sampling 
the non-local potential using the techniques of Fahy et 
al. |§. 

We have optimized wave functions using both the 
Ewald and MPC interactions. The wave functions gener- 
ated using the different interactions are almost identical. 
Properties other than the energy, such as pair correlation 
functions, are therefore hardly affected by the choice of 
interaction. As the MPC interaction gives the correct in- 
teraction between the electrons at short distances it may 
give a better account of, for example, the short distance 
behaviour of the pair correlation function, but more nu- 
merical work is required to investigate this point. 

In Fig. H we show results for VMC calculations of the 
energies of the same systems as in the previous section. 
These results are similar to those given in Ref. but 
they have been recalculated for this work with more ac- 
curate wave functions and better sampling, and the data 
have been corrected for the IPFSE. The total energies 
were calculated to a statistical accuracy of ±0.01 eV 
per atom. The VMC data display a smaller CFSE than 
the HF data, probably because the HF exchange hole 
is longer ranged than the screened hole obtained in the 
correlated calculations. For n = 2, the MPC interaction 
reduces the VMC finite size error by more than 50%, 
from -0.403 to -0.187 eV per atom. 

C. Application to DMC 

In the DMC method (|||, imaginary time evolution 
of the Schrodinger equation is used to evolve an ensem- 
ble of 3iV-dimensional electronic configurations towards 
the ground state. The calculations are made tractable 
by using the fixed node approximation and by incorpo- 
rating importance sampling. The method generates the 
distribution </>t^>) where tp is the best (lowest energy) 
wave function with the same nodes as the guiding wave 
function, <p^ . The accuracy of the fixed node approxima- 
tion can be tested on small systems and the results are 
normally very satisfactory. Q 

We evaluated the non-local pseudopotential energy us- 
ing the "locality approximation" . |2j| The short-time ap- 
proximation for the Green's function was used with a 
time step of 0.01 a.u. Li et al. |2(| found that using a 
time step of 0.015 a.u. gave a time-step error of less than 
0.03 eV per atom in silicon, so our time step error should 
be even smaller. The total energies were calculated to a 
statistical accuracy of ±0.02 eV per atom. 

Our MPC interaction is more complicated to apply in 
DMC than in VMC because the evaluation of the im- 
portance sampled Green's function requires the local en- 
ergy. The modified Hamiltonian, and hence the local 



energy, depends on the charge density, and therefore we 
must know the charge density before we can perform the 
DMC calculation. Fortunately, however, the local energy 
is relatively insensitive to the charge density used in the 
Hamiltonian (Eq. ||) because v^(r) — /(r) is small when 
|r| is small. 

We have tested the sensitivity of the Green's function 
to the charge density used in the Hamiltonian. Two 
candidate charge densities are the charge density of the 
determinantal part of the QMC wave function and the 
charge density of the VMC guiding wave function. Even 
for a small system (n = 2) we find it sufficient to use 
the LDA charge density during the calculation of the 
Green's function and to re-evaluate the charge density 
dependent term in the interaction energy using the DMC 
density obtained at the end of the calculation. The sen- 
sitivity rapidly reduces with increasing system size, and 
this procedure gives errors of less than 0.03 eV per atom 
for n — 2, and less than 0.01 eV per atom for larger sys- 
tem sizes. Therefore, the requirement of having a good 
approximation to the charge density in advance of the 
DMC calculation does not pose a significant difficulty. A 
successful DMC calculation requires a good quality VMC 
trial function and its charge density can be obtained dur- 
ing the process of wave function optimization. 

In Fig. | we show results for DMC calculations on the 
same systems as for our VMC study, the largest of which 
contains 1000 electrons. The results include a correction 
for the IPFSE. The convergence behaviour is very sim- 
ilar to the VMC data. The MPC energies are always 
above those for the Ewald interaction and the MPC in- 
teraction significantly reduces the CFSE. These results 
demonstrate that the finite size errors within VMC and 
DMC calculations are very similar and that the MPC 
interaction is similarly effective in both methods. 

We have fitted the residual CFSE errors in both the 
VMC and DMC calculations to the form b/N x . The fit is 
reasonable and gives a value of x close to unity for both 
the Ewald and MPC data. For these data it is possi- 
ble to obtain more accurate approximations to E^ AG by 
using such an extrapolation function, although the cal- 
culations for the large system sizes are costly, especially 
within DMC. The extrapolated energies should be more 
accurate for the MPC interaction because the extrapola- 
tion term is significantly smaller. 

Many interesting applications of VMC and DMC 
methods will be to problems in which the quantity of 
physical interest is the difference in energy between two 
large systems. Examples of such problems are calcula- 
tions of excitation energies and defect energies in solids. 
In such cases the energy of interest is approximately inde- 
pendent of the size of the simulation cell, so that for each 
simulation cell size it is the energy of the whole simulation 
cell which must be converged to the required tolerance, 
not the energy per atom as we plotted in Figs. ^ and ^. 
In these cases extrapolation will be so costly that it can 
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hardly be contemplated. In some cases the CFSE will 
largely cancel between the two systems, as occurs in our 
excitation energy calculations (see next section). This 
cancellation cannot always be relied upon, however, espe- 
cially when the simulation cells contain different numbers 
of particles, and the use of the MPC interaction should 
be particularly beneficial in such cases. 

VII. EXCITATION ENERGIES 

The quasiparticle excitation energies are the energies 
for adding an electron to the system or subtracting one 
from it. A quasiparticle energy has both real and imag- 
inary parts, the latter giving the quasiparticle lifetime. 
For the minimum energy electron and hole quasiparti- 
cle excitations the imaginary parts of the quasiparticle 
energies are zero and the quasiparticles have an infinite 
lifetime. In this case the exact quasiparticle energy gap 
can be written as 

Eg = (En+i — En) + (En -i — En) , (18) 

where En+i, En-i, and En are the ground state total 
energies of the N + 1, N — 1 and ./V electron systems. 
A general quasiparticle energy gap cannot be written in 
terms of differences between energies of exact eigenstates 
of the system, but such an approximation is often ac- 
curate for low energy gaps. The quasiparticle energies 
are measured in photoemission and inverse photoemis- 
sion experiments. In an optical absorption experiment 
a different process occurs in which an electron is excited 
from the valence to the conduction band. This introduces 
two quasiparticles into the system, the electron and hole, 
which interact and can form an exciton, in which case 
the lowest excitation energy is smaller than E g by the 
exciton binding energy, E\,. 

The HF method gives approximations to the energies 
of quasiparticles and the interactions between them. Ac- 
cording to Koopmans' theorem, if orbital relaxation is 
neglected, the quasiparticle energies are equal to the HF 
eigenvalues. Koopmans' theorem can be extended to 
include correlation effects. |27],|2^] The extended Koop- 
mans' theorem has been used in conjunction with VMC 
methods to calculate quasiparticle energies in silicon. |2^] 
In both of these methods the "quasiparticle energies" are 
real as they are obtained as approximations to the energy 
differences between exact eigenstates of the system. 

Recently there has been significant progress in apply- 
ing QMC techniques to calculate approximate excitation 
energies from eigenstates using "direct methods" . In 
these approaches an excitation energy is obtained by per- 
forming separate QMC calculations for the ground and 
excited states. A Slater- Jastrow wave function is used 
for the ground state, and for the optical gap the ex- 
cited state is formed by replacing a valence band single 



particle orbital by a conduction band one. We call this 
a "promotion" calculation, and such calculations have 
been reported for a nitrogen solid pOfl , diamond pl|p2[ , 
and silicon p3fl . Photoemission/inverse photoemission 
gaps may be obtained by using QMC to calculate the 
ground state energies of the N + 1, N — 1, and ./V elec- 
tron systems. Wave functions for the N + 1 and N — 1 
electron systems may be formed by adding or subtract- 
ing an orbital from the up- or down-spin determinants 
of the TV-electron wave function. We call this an "addi- 
tion/subtraction" calculation. For calculations with pe- 
riodic boundary conditions the simulation cell is made 
charge neutral by adding a compensating uniform back- 
ground charge density. Calculations of this type have 
been reported for one- [Q and two-dimensional JTJ] 
model systems, while results for a three-dimensional sys- 
tem (silicon) are reported in this paper. 

QMC calculations of excitation energies in extended 
systems are computationally very demanding because 
they are '-^' effects, i.e., the fractional change in energy 
is inversely proportional to the number of electrons in 
the system. This means that high statistical accuracy 
is required to obtain good results. The largest system 
for which excitation energies have been calculated prior 
to this paper is 16 atoms (64 electrons). [[33| The total 
finite size error in the ground state energy for that sys- 
tem was estimated to be about 16 eV per simulation cell, 
while the energy scale of interest for the excitations is 
of order 0.1 eV. Like almost all methods for calculating 
excitation energies, QMC calculations of this type only 
work because of a strong cancellation of errors between 
the ground and excited states. It turns out that the finite 
size errors tend to reduce the energy gap, while the errors 
in the trial wave functions are usually larger for the ex- 
cited states than for the ground state and so increase the 
energy gap. Although good agreement with experimental 
excitation energies has been found using small simulation 
cells, [|3l|-|33}| one is left with the suspicion that if larger 
simulation cells were used the agreement with experi- 
ment might be significantly worse because the finite size 
effects would be smaller. Before these QMC techniques 
can be relied upon for calculating excitation energies it is 
necessary that the issue of finite size effects be properly 
explored. In the next sections we address the following 
questions: 

1. What are the sizes and origins of the finite size 
effects in excitation energy calculations? 

2. What are the differences in finite size effects be- 
tween promotion and addition/subtraction calcu- 
lations? 
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A. HF theory of excitation energies 

First we consider excitation energies for solids within 
HF theory. The HF equations with the MPC interaction 
are obtained by minimizing the HF energy of Eq. [L6| with 
respect to the single-particle orbitals, giving 

-\v 2 4>,+ J p(r> E (r-r')dr' 
- E S «>sJ W^(r')/(r - r') dr' 



%3 

+ Vext<j>i = U 



(19) 



If we neglect the relaxation of the orbitals, the energy 
required to excite an electron from the j th (occupied) 
orbital into the i th (unoccupied) orbital is 



AEij = (a - ej) 



J pi(r)pj (t')ve(t - r')dr dr' 

+ S SzSj J ^(r)0*(r')^-(r)^(r')/(r - r')drdr' 

+ i f p i (r)p i (r')[v E (r-r')-f(r-r')} drdr' 

+ \j Pj(')PiP) iMr - r') - /(r - r')] drdr' , (20) 

where pk — \4>k\ 2 is the charge density from the k th 
orbital. The first term is the eigenvalue difference for 
the excitation while the second and third terms are the 
Hartree and exchange interactions between the electron 
and hole. Within this approximation the electron-hole 
terms go to zero in the limit of an infinitely large simu- 
lation cell. The fourth and fifth terms on the right hand 
side are absent if one uses the Ewald interaction instead 
of the MPC interaction, i.e., we replace / by i>e- When 
the relaxation of the orbitals is neglected these terms also 
go to zero when the size of the simulation cell goes to in- 
finity because we tends to 1 fr over most of the simulation 
cell. 

The addition/subtraction gap is given by 
E g = (Elf E^) (£„ HF - E™) 

+ \J Pi(r)Pi(r') Mr - r') - /(r - r')] drdr' 

+ \f Pi(r) Pj (T f ) Mr - r') - /(r - r')] drdr' , (21) 

where E^ F is the HF ground state energy of the N- 
electron system, E^J is the energy of the state with an 
electron added to the i th (previously unoccupied) orbital, 
along with the uniform background charge, and E_J is 
the energy of the state where an electron is removed from 
the j th orbital, along with the background charge. The 



standard Koopmans' theorem has been modified and con- 
tains two additional terms, which also occur in the pro- 
motion energy, AEij . We have evaluated these additional 
terms using LDA orbitals and have found that even for a 
small simulation cell (n = 2) they are very small, being 
in the range ±0.05 eV, and they decrease rapidly with 
system size. We do not expect that the use of exact HF 
orbitals or orbital relaxation will greatly affect these re- 
sults. 

In Fig. H we show the addition/subtraction energies 
calculated using Eq. ^l] for the T25' — > T15 energy gap 
and the valence band width, calculated with both the 
Ewald and MPC interactions, along with LDA values. 
(The results for other energies show similar behaviour.) 
We do not show the promotion energies in Fig. | be- 
cause they differ from the addition/subtraction energies 
only by the exciton binding energy, which decreases with 
increasing system size quite rapidly. Fig. || shows that 
the HF results for the Ewald and MPC interactions are 
very similar. The band width converges by about n = 
7, but the band gap is still slowly increasing at n = 12, 
and the Ewald and MPC values are not yet equal, which 
they must be at convergence. For the largest system 
size studied (n = 12) the MPC gap and valence band 
width are 7.4 eV and 17.7 eV respectively, which are 
a little smaller than the HF values of 8.0 eV and 18.9 
eV given in Ref. |55|. Presumably the major reasons for 
these differences are that we use LDA wave functions and 
LDA-derived pseudopotentials, although as noted above 
there is clear evidence that in our calculations the HF 
energy gap has still not fully converged at n = 12. The 
LDA excitation energies converge very rapidly with sys- 
tem size. Note that this would not be true in either LDA 
or HF theory if we studied isolated clusters of atoms. 
In a recent study of silicon clusters, Ogiit et al. ||(J 
found large differences between the band gap in the LDA 
eigenvalues and the band gap calculated by electron addi- 
tion/subtraction. As shown by Franceschetti et al. [p7| , 
these differences are due to the charging of the cluster 
when an electron is added or subtracted, which does not 
occur in our calculations because a uniform background 
is added to preserve charge neutrality. The slow conver- 
gence of the HF excitation energies apparent in Fig. || 
therefore arises from the exchange energy. Moreover, be- 
cause the results with the Ewald and MPC interactions 
are almost the same, the source of the error is not the 
interaction with the exchange hole, but the shape of the 
exchange hole. This is because the excitation energy de- 
pends on the change in the exchange hole due to the 
excitation, which is not strongly localised. 

In a very interesting set of calculations Engel et al. |l4| ] 
studied excited states of a model two-dimensional system 
using LDA, GW, VMC and DMC techniques. They per- 
formed a number of VMC calculations with increasing 
system size and found that the addition/subtraction gap 
tended to increase with system size, which is the same 
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behaviour as we have found in our HF calculations. En- 
gel et al. went on to give an explanation of this effect. 
Their explanation was that in the N + 1 (N — 1) electron 
systems there is an additional electrostatic energy due 
to the interaction of the extra electron (hole) with the 
additional electron (hole) in the other simulation cells. 
Taking into account the additional compensating uni- 
form background charge that was added to keep each cell 
neutral, this additional energy is negative and therefore 
the energies £jv+i and £jv_i are lower than they should 
be. Engel et al. showed that the observed finite size ef- 
fect is much smaller than the Madelung energy for point 
charges, and to explain this they argued that the effect 
would be screened by the response of the other electrons. 
This argument implies that the finite size effects in ad- 
dition/subtraction calculations are larger than those in 
promotion calculations. 

Our analysis of the situation is as follows. For simplic- 
ity we consider our HF calculations, where the interaction 
energy can be divided into Hartree and exchange contri- 
butions. The significant underestimation of the HF band 
gaps of small systems is not due to the Hartree terms, 
which by construction are the same as for our LDA cal- 
culations and give very small finite size effects in the band 
gaps. The finite size error in the HF gaps therefore arises 
from the exchange energy. By comparing band gaps cal- 
culated with the Ewald and MPC interactions we can see 
whether the problem lies with the interaction or with the 
shape of the exchange hole. Because we find that the 
band gaps calculated with the Ewald and MPC interac- 
tions are very similar we conclude that the form of the 
interaction is not the important consideration. Therefore 
the source of the problem must be the finite size errors 
present in the shape of the exchange hole. This argument 
implies that the finite size effects in addition/subtraction 
calculations are similar to those in promotion calcula- 
tions. Our viewpoint is supported by the HF results that 
have been presented in this sub-section and also by the 
VMC results to be discussed in the next sub-section. 

In summary, the HF excitation energies calculated 
with the Ewald and MPC interactions are very similar. 
Within HF theory the largest finite size error in excita- 
tion energies arises from the shape of the exchange hole, 
which leads to slow convergence with system size. The 
finite size errors in promotion and addition/subtraction 
HF calculations are of similar size. 

B. QMC theory of excitation energies 

We now apply the theory developed in the previous 
section to QMC calculations of excitation energies. Al- 
though we have just demonstrated that the HF gaps con- 
verge rather slowly with system size, we showed earlier 
that the finite size effects in VMC and DMC ground state 



energies are smaller than in HF theory. It is important 
to see whether this also applies to excitation energies. 

VMC is computationally cheaper than DMC and so we 
are able to compute excitation energies using VMC over 
a larger range of system sizes. We expect that the finite 
size effects in DMC will follow those in VMC, as our VMC 
calculations retrieve about 90% of the fixed-node correla- 
tion energy. We have computed the T25' ^Ti5 excitation 
energy in silicon within VMC for the system sizes n — 
1, 2, 3, 4 using both promotion and addition/subtraction 
methods. The calculations were performed with k s = T 
and the other computational details were the same as for 
the ground state calculations. We used Jastrow factors 
optimized for the ground state of each system which were 
left unchanged for the excited state. In tests on the n = 
2 system we found that separately optimizing the Jas- 
trow factors for both ground and excited states did not 
significantly change the results. The computational cost 
of the n — 4 calculations is very large; an error bar in 
the excitation energies of ±0.3 eV requires an error bar 
of ±0.0006 eV per electron. Although the computational 
effort is large we believe that such a study is necessary 
to establish the accuracy of QMC excitation energy cal- 
culations. 

In Fig. U we show the excitation energies obtained 
with the Ewald interaction via the promotion and ad- 
dition/subtraction methods. (Results for the MPC in- 
teraction are very similar.) The promotion and addi- 
tion/subtraction results are nearly the same, but the pro- 
motion energies are slightly smaller because they include 
an exciton binding energy, which decreases as the sys- 
tem size increases. The results are consistent with a slow 
increase in the excitation energy with system size and in- 
dicate that reasonable convergence is already obtained at 
n — 2. The increase in excitation energies with system 
size is the same trend as in HF calculations, although 
the finite size errors are smaller in the correlated cal- 
culations. The finite size errors in the promotion and 
addition/subtraction methods are not significantly dif- 
ferent at this level of statistical accuracy. On general 
grounds we expect the finite size effect in the promotion 
calculation to be slightly larger. This follows because 
the trend for both promotion and addition/subtraction 
calculations is that the excitation energy is reduced for 
small system sizes and this effect is enhanced in the pro- 
motion calculation by the exciton binding energy which 
is larger for small systems. The calculated excitation 
energy is roughly 4 eV, which is larger than the experi- 
mental value of 3.40 eV |58|, and also a little larger than 
our DMC value for an n = 2 simulation cell of 3.7 eV |33| . 
Our study demonstrates that the largest contribution to 
the error in the VMC band gap for n > 2 arises from the 
approximate nature of the trial wave functions, and not 
from finite size effects. 

The exciton binding energy can be calculated 
as the difference between the promotion and addi- 
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tion/subtraction gaps. The exciton binding energy for 
the T25' — > Ti5 excitation is small and we can only re- 
solve it from the statistical noise for the smallest (n = 1) 
cell, which gives a value of 0.28 ±0.01 eV. In earlier QMC 
calculations |5(]-|32| the exciton binding energy was esti- 
mated using the Mott-Wannier formula 



where e is the relative permittivity and r is the radius of 
the localisation region. Using e = 11.7 and the appropri- 
ate radius for n = 1 of r = 4.0 a.u. gives an exciton bind- 
ing energy of 0.29 eV. This is extremely close to the VMC 
value, but the excellent agreement is probably fortuitous 
since the n = 1 cell is so small that it is appropriate to 
use a value of e at finite wave vector, which would be 
smaller. The exciton binding energy may also be eval- 
uated within HF theory as the sum of the second and 
third terms on the r.h.s. of Eq. [2(]. This gives 0.75 eV 
for the n = 1 cell using the Ewald interaction, which is 
considerably larger than the VMC value because the lat- 
ter calculation includes screening of the electron-electron 
interaction. 

Note that the promotion and addition/subtraction 
methods differ significantly in the required computa- 
tional effort. Suppose, for example, that we wish to 
calculate an energy gap by either the promotion or ad- 
dition/subtraction methods. Let us assume that the in- 
trinsic variance of the local energy is the same for each 
of the energies, which is a good approximation for our 
silicon calculations, and suppose that an acceptable er- 
ror bar is obtained in a promotion calculation by per- 
forming M Monte Carlo moves for both the ground and 
excited states. A simple calculation shows that the most 
efficient way to achieve the same error bar in an addi- 
tion/subtraction gap is to perform 2M moves for each of 
the N + 1 and N — 1 systems and 4M for the ground 
state, giving a total cost of 8M moves. It is therefore 
four times more expensive to calculate an energy gap to 
some given accuracy by the addition/subtraction method 
than by the promotion method. 

C. Modified interaction for excitation energies 

In this section we introduce a modified electron- 
electron interaction specifically designed to describe exci- 
tation energies within periodic boundary conditions sim- 
ulations. Two problems arise when trying to model 
excitations using finite simulation cells subject to peri- 
odic boundary conditions. One is that the excitation is 
"squeezed" into the simulation cell, and the other is that 
there are spurious interaction betweens the periodic repli- 
cas of the simulation cell. Here we address the problem 
of the spurious interactions using the ideas of the MPC 



interaction. With our MPC interaction the replicas in- 
teract only via the Hartree energy. The charge density 
on promotion or addition/subtraction of an electron can 
be written as the sum of the ground state charge den- 
sity and a small deviation, i.e., p(r) = p(r) + A(r). We 
can modify the Hartree term so that this charge density 
interacts with the ground state density in the replicas. 
This leads to the interaction energy 

J i>j k 

+ j p{v)p{t') [t> E (r - r') - /(r - r')] drdr' 

~ \ / P(r)p(r') K(r - r') - /(r - r')] drdr' , (23) 

where \4>\ 2 generates the charge density p and the ground 
state charge density, p, is fixed. A HF analysis of this 
interaction shows that the HF equations are identical 
to Eq. [n], so that the orbitals and eigenvalues are un- 
altered. However, the ground and excited state energy 
expressions are modified. For the excited states we ob- 
tain analogues of Eqs. [|(] and El], but without the terms 
involving (we — /), i.e., we retrieve the standard Koop- 
mans' theorem. We have already shown that these terms 
are small for silicon, although they will be significant 
in cases when the change in the charge density due to 
the excitation is strongly localized. This analysis pro- 
vides further evidence that the electrostatic interactions 
between the simulation cell and its replicas is not neces- 
sarily an important source of finite size error in excited 
state energy calculations. 

VIII. CONCLUSIONS 

Large Coulomb finite size errors arise in total energy 
calculations when using the Ewald form of the electron- 
electron interaction. These finite size errors may be 
greatly reduced by using our MPC interaction in which 
the inter-particle interaction is exactly equal to 1/r at 
short distances and the long range interactions are re- 
placed by a mean-field-like one-electron potential. It 
is consistent to use the MPC interaction in conjunction 
with "independent particle finite size corrections" derived 
from density functional calculations, as long as the latter 
calculations are performed with the exchange-correlation 
functional appropriate to the infinite system. Although 
the long range mean-field-like contribution to the MPC 
interaction involves the charge density of the system, to- 
tal energies are insensitive to its form and only an ap- 
proximate charge density is required. 

The MPC interaction can be used consistently for all 
Coulomb interactions in the system, although it normally 
reduces to using the MPC interaction for the electron- 
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electron interaction, retaining the standard Ewald inter- 
action for the electron-ion and ion-ion terms. The Ewald 
and MPC interactions may be used in tandem as an ef- 
ficient diagnostic of Coulomb finite size errors. If the 
Ewald and MPC results agree then the Coulomb finite 
size error should be small. 

If the simulation cell is too small then the confinement 
of the XC hole makes the XC energy more negative. This 
source of error is intrinsic to using a finite simulation cell. 
However, even when the XC hole is artificially confined 
by a small simulation cell the MPC interaction still gives 
a better estimate of the XC energy than the Ewald in- 
teraction. 

Excitation energies calculated within fixed-LDA- 
orbital HF theory show significant finite size effects. 
However, in correlated calculations the finite size ef- 
fects are smaller and accurate excitation energies can 
be obtained using quite small simulation cells. In sili- 
con we find that the finite size errors in VMC electron- 
promotion ("optical absorption") and electron addi- 
tion/subtraction ( "photoemission" ) calculations are sim- 
ilar, and that the optical promotion method has the 
greater statistical efficiency. The finite size errors for low 
lying excitations in silicon are small, and quite accurate 
results may be obtained from 16 atom cells. 

We have described new developments aimed at un- 
derstanding and reducing finite size errors in many- 
body quantum simulations using periodic boundary con- 
ditions. Since one cannot get exact answers for an infi- 
nite system from a finite simulation cell whatever inter- 
action is used, there is no "exact interaction" for a finite 
system with periodic boundary conditions. The Ewald 
and MPC interactions are alternative model interactions 
compatible with periodic boundary conditions, and the 
relevant question is which model interaction gives results 
which most closely approximate those for very large sim- 
ulation cells. The Ewald and MPC interactions differ 
by an amount which is inversely proportional to the size 
of the simulation cell and therefore they give the same 
energy per particle in the limit of an infinitely large simu- 
lation cell. However, for finite cells the Ewald and MPC 
interactions can give significantly different energies. In 
every test we have performed the energy calculated with 
the MPC interaction is closer than the Ewald energy to 
the value for a very large system. The MPC interaction 
is applicable to both metals and insulators and it is faster 
to compute than the Ewald interaction. Given these facts 
we believe that the MPC interaction should be used for 
all quantum many-body calculations of total energies in 
systems with periodic boundary conditions. 
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FIG. 1. Convergence of the LDA and HF ground state en- 
ergies per atom of silicon as a function of simulation cell size, 
n, using F-point BZ sampling. 

FIG. 2. The LDA and HF ground state energies per atom 
of silicon as a function of simulation cell size, n, using L-point 
BZ sampling. 

FIG. 3. The VMC ground state energies per atom of silicon 
as a function of simulation cell size, n. A correction for the 
IPFSE is included. The statistical error bars are ±0.01 eV 
per atom. 

FIG. 4. The DMC ground state energies per atom of silicon 
as a function of simulation cell size, n. A correction for the 
IPFSE is included. The statistical error bars are ±0.02 eV 
per atom. 

FIG. 5. Convergence of the T 2 5' — * Ti5 excitation energy 
and valence band width of silicon. The data correspond to 
addition/subtraction energies calculated within HF theory as 
a function of simulation cell size, n, using both the MPC and 
Ewald interactions. LDA results are also shown. 

FIG. 6. The r 25 ' — > Fis excitation energy of silicon calcu- 
lated within VMC theory as a function of simulation cell size, 
n. Data for the Ewald interaction are shown obtained via 
both the promotion and addition/subtraction methods. 
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Figure 1 Kent et. al. 




Figure 2 Kent et. al. 
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Figure 3 Kent et. al. 
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Figure 4 Kent et. al. 
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